function D_matrixs(E::Float64, μ::Float64)
	D0 = E/(1.0 - μ*μ)
	return D0*[1.0 μ 0.0; μ 1.0 0.0; 0.0 0.0 0.5*(1.0 - μ)]
end

function D_matrixc(E::Float64, G::Float64)
	return [E 0.0; 0.0 5.0*G/6.0]
end

function ksteel(xi::Float64, xj::Float64)
   k = zeros(Float64, ndof, ndof)
   det = (xj - xi)/2.0

   a1 = -0.5*hweb - fbot; b1 = -0.5*hweb;
   a2 = -0.5*hweb; b2 = 0.5*hweb;
   a3 = 0.5*hweb; b3 = 0.5*hweb + ftop;
   D0 = D_matrixs(Es, μs)

   # I-shaped steel beam integration
   # part 1, lower flange
   for i = 1:ngpx # integrate in x axis
        for j = 1:ngpz # integrate in z axis
            z = 0.5*(b1 - a1)*gaussdataz[j, 1] + 0.5*(b1 + a1)
            w = gaussdatax[i, 2]*0.5*(b1 - a1)*gaussdataz[j, 2]
			B = BMatrixS(xi, xj, i, z)
			int = (B') * D0 * B
            k = k + int*det*wbot*w
        end
    end
	# part 2, web
	for i = 1:ngpx # integrate in x axis
		 for j = 1:ngpz # integrate in z axis
			 z = 0.5*(b2 - a2)*gaussdataz[j, 1] + 0.5*(b2 + a2)
			 w = gaussdatax[i, 2]*0.5*(b2 - a2)*gaussdataz[j, 2]
			 B = BMatrixS(xi, xj, i, z)
			 int = (B') * D0 * B
			 k = k + int*det*wweb*w
		 end
	 end
	 # part 3, upper flange
	 for i = 1:ngpx # integrate in x axis
		  for j = 1:ngpz # integrate in z axis
			  z = 0.5*(b3 - a3)*gaussdataz[j, 1] + 0.5*(b3 + a3)
			  w = gaussdatax[i, 2]*0.5*(b3 - a3)*gaussdataz[j, 2]
			  B = BMatrixS(xi, xj, i, z)
			  int = (B') * D0 * B
			  k = k + int*det*wtop*w
		  end
	  end
	return k
end

function kcs(xi::Float64, xj::Float64)
	k = zeros(ndof,ndof)
	det = (xj - xi)/2.0

	for i = 1:ngpx
			w = gaussdatax[i, 2]
			n1 = Nucs(i, h1, h2)
			n2 = Nwcs(i, h2)
			int = n1' * kucs * n1 + n2' * kwcs * n2
			k = k + int*det*w
	end
	return k
end

function krebar(xi::Float64, xj::Float64)
	# z: rebar possition in o-zs coordinate
	k = zeros(ndof, ndof)
	ngpx = size(gaussdatax,1)
	det = (xj - xi)/2.0
	for i = 1:ngpx
		w = gaussdatax[i, 2]
		n = Nϵc(xi, xj, i, zr)
		int = n' * n
		k = k + int*Ar*(Er-Ec)*det*w
	end
	return k
end

function kconcrete(xi::Float64, xj::Float64)
	k = zeros(ndof, ndof)
	det = (xj - xi)/2.0
	D0 = D_matrixc(Ec, Gc)
	for i = 1:ngpx
		for j = 1:ngpz
			z = 0.5*Hc*gaussdataz[j, 1]
			w = gaussdatax[i, 2]*0.5*Hc*gaussdataz[j, 2]
			B = BMatrixC(xi, xj, i, z)
			int = (B') * D0 * B
			k = k + int*det*Bc*w
		end
	end
	return k
end

function R0(xi::Float64, xj::Float64)
	det = (xj - xi)/2.0;
	r = zeros(ndof)
	for j = 1:ngpx
		w = gaussdatax[j, 2]
		r = r + det*w*reshape(Nwc(j, 0), ndof)*q0
	end
	return r
end

function R1(nel, nt)
	int = zeros(ndof);
	num = gnum[nel,:]
	elemcoord = gcoord[num,:]
	xi = elemcoord[1]; xj = elemcoord[3]
	det = (xj - xi)/2.0
	for i = 1:ngpz
		for j = 1:ngpx
			z = 0.5*Hc*gaussdataz[i, 1]
			w = gaussdatax[j, 2]*0.5*Hc*gaussdataz[i, 2]
			s = 0.0
			for i = 1:(nt-1)
				s = s + psi(nt, i, time0)*fem_dot(StressHistory[i, nel, :], Nstress(j, z))
			end
			Ks = reshape(Nϵc(xi, xj, j, z),ndof)*(en(nt, time0)*1e6*epsilon_shn(nt, time0) - s)
			int = int + Bc*det*w*Ks
		end
	end
	for j = 1:ngpx
		w = gaussdatax[j, 2]
		s = 0.0
		for i = 1:(nt-1)
			s = s + psi(nt, i, time0)*fem_dot(StressHistory[i, nel, :], Nstress(j, zr))
		end
		Ks = reshape(Nϵc(xi, xj, j, zr),ndof)*(s - en(nt, time0)*1e6*epsilon_shn(nt, time0))
		int = int + Ar*det*w*Ks
	end

for i=1:ngpz
	for j = 1:ngpx
		z = 0.5*Hc*gaussdataz[i, 1]
		w = gaussdatax[j, 2]*0.5*Hc*gaussdataz[i, 2]
		s = 0.0
		for i = 1:(nt-1)
			s = s + psi(nt, i, time0)*fem_dot(StressHistorytao[i, nel, :], Nstress(j, zr))
		end
		Ks = -reshape(Nγc(xi, xj, j), ndof)*s
		int = int + Bc*det*w*Ks
	end
end

	return int
end

function Nstress(ξi, z)
	nx = zeros(ndp); nx[ξi] = 1.0;
	nz = lagrange(0.5*Hc*gz, z)
	np = zeros(3*ndp)
	count = 0
	for iz = 1:3
		for ix = 1:ndp
			count = count + 1
			np[count] = nx[ix]*nz[iz]
		end
	end
	return np
end

function strainc(nel, sol, ξi, z)
   elemcoord = gcoord[gnum[nel,:],:]
   xi = elemcoord[1]
   xj = elemcoord[ndp]
   de = [extract(sol, g_g[k, nel]) for k = 1:ndof]
   return fem_dot(Nϵc(xi, xj, ξi, z), de)
end

function strainctao(nel, sol, ξi, z)
   elemcoord = gcoord[gnum[nel,:],:]
   xi = elemcoord[1]
   xj = elemcoord[ndp]
   de = [extract(sol, g_g[k, nel]) for k = 1:ndof]
   return fem_dot(Nγc(xi, xj, ξi), de)
end

function stressc(sol, nel, ξi, z, nt)
	#(* sol is the solution of the ntth time tick *)
	s = 0.0
	for i = 1:(nt-1)
		s = s + psi(nt, i, time0)*fem_dot(StressHistory[i, nel,:], Nstress(ξi, z))
	end
	r = en(nt, time0)*1e6*(strainc(nel, sol, ξi, z) - epsilon_shn(nt, time0)) + s
	return r
end

function stressctao(sol, nel, ξi, z, nt)
	#(* sol is the solution of the ntth time tick *)
	s = 0.0
	for i = 1:(nt-1)
		s = s + psi(nt, i, time0)*fem_dot(StressHistorytao[i, nel,:], Nstress(ξi, z))
	end
	gc = (5.0/6.0) * en(nt, time0)/(2.0*(1 + μc)) * 1e6
	r = gc*strainctao(nel, sol, ξi, z) + s
	return r
end

function record(sol, nt)
	#   (* sol is the solution of the ntth time tic *)
	sel = zeros(nels, 3*ndp)
	seltao = zeros(nels, 3*ndp)
	for nel = 1:nels
		count = 0
		for iz = 1:3
			for ix = 1:ndp
				count = count + 1
				z = 0.5*Hc*gz[iz]
				sel[nel, count] = stressc(sol, nel, ix, z, nt)
				seltao[nel, count] = stressctao(sol, nel, ix, z, nt)
			end
		end
	end
	StressHistory[nt,:,:] = sel
	StressHistorytao[nt,:,:] = seltao
	return nothing
end

function assemble0(gnum::Array{Int64,2}, gcoord::Array{Float64,2}, nf::Array{Int64,2}, nels::Int64, neq::Int64)
    gK = zeros(Float64, neq, neq)
    gR = zeros(Float64, neq)

    for iel = 1:nels
        num = gnum[iel,:]
        elemcoord = gcoord[num,:]
		xi = elemcoord[1]; xj=elemcoord[ndp];
        g = numtog(num, nf)
        kma = ksteel(xi, xj) + kcs(xi, xj)
		rvec = R0(xi, xj)
	    for i = 1:length(g)
			row = g[i]
			if row != 0
				gR[row] = gR[row] + rvec[i]
			end
	        for j = 1:length(g)
	            col = g[j]
	            if row*col != 0
	                gK[row, col] = gK[row, col] + kma[i, j]
	            end
	        end
	    end
	end

	return (gK, gR)
end

function assemble1(gnum::Array{Int64,2}, gcoord::Array{Float64,2}, nf::Array{Int64,2}, nels::Int64, neq::Int64)
    gK = zeros(Float64, neq, neq)
    gR = zeros(Float64, neq)

    for iel = 1:nels
        num = gnum[iel,:]
        elemcoord = gcoord[num,:]
		xi = elemcoord[1]; xj=elemcoord[ndp];
        g = numtog(num, nf)
        kma = krebar(xi, xj) + kconcrete(xi, xj)
		rvec = R1(iel, nt)
	    for i = 1:length(g)
			row = g[i]
			if row != 0
				gR[row] = gR[row] + rvec[i]
			end
	        for j = 1:length(g)
	            col = g[j]
	            if row*col != 0
	                gK[row, col] = gK[row, col] + kma[i, j]
	            end
	        end
	    end
	end

	return (gK, gR)
end

function assemble_loads0(gnum::Array{Int64,2}, gcoord::Array{Float64,2}, nf::Array{Int64,2}, nels::Int64, neq::Int64)
    gR = zeros(Float64, neq)
    for iel = 1:nels
        num = gnum[iel,:]
        elemcoord = gcoord[num,:]
		xi = elemcoord[1]; xj=elemcoord[ndp]
        g = numtog(num, nf)
		rvec = R0(xi, xj)
	    for i = 1:length(g)
			row = g[i]
			if row != 0
				gR[row] = gR[row] + rvec[i]
			end
	    end
	end

	return gR
end

function assemble_loads1(gnum::Array{Int64,2}, gcoord::Array{Float64,2}, nf::Array{Int64,2}, nels::Int64, neq::Int64)
    gR = zeros(Float64, neq)
    for iel = 1:nels
        num = gnum[iel,:]
        g = numtog(num, nf)
		rvec = R1(iel, nt)
	    for i = 1:length(g)
			row = g[i]
			if row != 0
				gR[row] = gR[row] + rvec[i]
			end
	    end
	end
	return gR
end
